Hydrodynamic interactions of spherical particles in Poiseuille flow between two 

parallel walls 



S. Bhattacharya, 1 J. Blawzdziewicz, 1 and E. Wajnryb 1,2 

department of Mechanical Engineering, Yale University, P.O. Box 20-8286, New Haven, CT 06520 

On leave from IPPT Warsaw, Poland 
(Dated: February 2, 2008) 

We study hydrodynamic interactions of spherical particles in incident Poiseuille flow in a channel 
with infinite planar walls. The particles are suspended in a Newtonian fluid, and creeping- flow condi- 
tions are assumed. Numerical results, obtained using our highly accurate Cartesian-representation 
algorithm [Physica A xxx, xx, 2005], are presented for a single sphere, two spheres, and arrays 
of many spheres. We consider the motion of freely suspended particles as well as the forces and 
torques acting on particles adsorbed at a wall. We find that the pair hydrodynamic interactions 
in this wall-bounded system have a complex dependence on the lateral interparticle distance due 
to the combined effects of the dissipation in the gap between the particle surfaces and the back- 
flow associated with the presence of the walls. For immobile particle pairs we have examined the 
crossover between several far-field asymptotic regimes corresponding to different relations between 
the particle separation and the distances of the particles from the walls. We have also shown that 
the cumulative effect of the far-field flow substantially influences the force distribution in arrays of 
immobile spheres. Therefore, the far-field contributions must be included in any reliable algorithm 
for evaluating many-particle hydrodynamic interactions in the parallel- wall geometry. 



I. INTRODUCTION 

In his pioneering work (more than eighty years ago) 
Faxen pj considered motion of a spherical particle sus- 
pended in a fluid confined by two parallel walls. A re- 
cent, considerable interest in particle motion in confined 
geometries has been stimulated by development of new 
experimental techniques @i S El S El an d by emerg- 
ing applications, such as the microfluidic devices and 
technologies for production of microstructured materials 
by a self-assembly process jj| EJj . 

There have been published a number of fundamental 
experimental and numerical studies on particle dynamics 
in channels with par allel planar walls for suspensions of 
Brownian El ETllBElJll El E EEcj and non- 
Brownian |2ll l22t l23t l24l l25j particles. Some of these 
studies focused on quasi-two-dimensional phenomena 0, 
El El El El Il9| , and some on bulk properties, such as 
particle migration in the pressure-driven [2l], [22, [23, |25j 
or shear j24j flow. 

Quantitative numerical studies of wall-bounded sus- 
pensions require efficient methods for evaluation of multi- 
particle hydrodynamic interactions in these systems. 
Some interesting numerical results were obtained with 
the help of the wall-superposition approximation [Tl Il4j 
or by m odelin g^ th e walls as arrays of immobile spheres 
[Il^lMllllll. These approaches seem sufficient for 
describing certain qualitative features of wall-bounded 
suspensions (e.g., in Stokesian-dynamic simulations of 
hydrodynamic particle diffusion) but the accuracy of 
such approximations is often unknown. Moreover, in 
some cases, they entirely miss certain important phe- 
nomena. For example, the superposition approximation 
reproduces neither the large transverse resistance coef- 
ficient of rigid arrays of spheres j2|J nor the enhanced 
relative transverse particle motion, observed by Cui et 



al. 1 191 and independently predicted by our recent anal- 
ysis [27]. 

To overcome these difficulties, we have developed an 
accurate Cartesian-representation method for evalua- 
tion of multiparticle hydrodynamic interactions in wall- 
bounded suspensions of spheres [2|| . (A related approach 
was also independently proposed by Jones [28].) Our 
method relies on expanding the flow in a wall-bounded 
system using two basis sets of Stokes flows. The spherical 
basis set of multipolar flows is used to describe the inter- 
action of the fluid with the particles, and the Cartesian 
basis set is used to account for the presence of the walls. 

In our previous studies, the Cartesian-representation 
method was applied to determine the resistance func- 
tions for systems of spheres in quiescent fluid [2II2I |29| . 
In the present paper we extend these results to suspen- 
sions in a pressure-driven external flow. We note that the 
one-particle motion in such a system was investigated by 
Jones and Staben et al. 30] but, to our knowledge, 
no accurate multiparticle results have been reported so 
far. 

This paper is organized as follows. In Sec.[n]the system 
is defined, and in Sec. 1 1 i 1 1 the Cartesian-representation 
method is summarized. Our numerical results for single- 
particle, two-particle, and multiparticle systems are de- 
scribed in Sec. II VI Conclusions are presented in Sec. [VJ 



II. PARTICLES IN PARABOLIC FLOW 

We consider a suspension of N spherical particles of 
diameter d = 2a in creeping flow between two parallel 
planar walls. The no-slip boundary conditions are sat- 
isfied on the walls and the particle surfaces. The walls 
are in the planes z = and z = H, where H denotes 
the wall separation, and (x, y, z) are the Cartesian coor- 
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dinates. The position of the center of particle i (where 
i = 1, . . . , N) is denoted by Rj, and its translational and 
rotational velocities are Uj and Slj, respectively. The ex- 
ternal forces and torques acting on particle i are denoted 
by Ti and Tj. 

In this paper we focus on particle motion in an imposed 
parabolic flow of the form 




where U p is the flow amplitude, and e x is the unit vec- 
tor along the x coordinate. The forces and torques on 
immobile particles with 

U, =0, Hi = 0, (2) 

can be represented by the resistance formula 
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Similarly, the velocities of freely suspended particles with 



the translational-rotational mobility matrix fi for a sys- 
tem of spherical particles confined between two paral- 
lel walls. In the present paper our method is used to 
evaluate the friction and mobility matrices <^f p and vf p 
(A = t, r) associated with the Poiseuille flow between the 
walls. 



III. CARTESIAN AND HELE-SHAW 
REPRESENTATION METHODS 

In this section we summarize the key elements of our 
Cartesian-representation method for evaluating the hy- 
drodynamic friction and mobility matrices in a suspen- 
sion confined between two parallel walls. We also outline 
our asymptotic results, which rely on expansion of the far 
field flow into a Hele-Shaw basis. The asymptotic results 
apply for sufficiently large interparticle separations. 

A detailed description of our technique is presented in 
|26j and in j29j. The Hele-Shaw basis and its relation 
to the spherical basis j^] used in our analysis [2(| is 
summarized in Appendices 1X1 and IB1 



Ti = 0, Ti = 0, 
can be represented by the mobility formula 
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In our considerations we assume that the applied flow 
JU is in the x direction. Thus, the resistance coefficients 
C- p and Ci P an d the mobility coefficients ia p and i/ p are 
vectors. For the external parabolic flow applied in an 
arbitrary lateral direction, the corresponding resistance 
and mobility coefficients have a tensorial character. 

Condition @ can be obtained by applying to immo- 
bile particles the forces and torques opposite to those 
given by equation ©. Thus, the resistance and mobility 
coefficients C and v satisfy the following relation 



N 

E 



(6) 



where u^ B (A,B = t,r) are the translational and ro- 
tational components of the usual mobility matrix |31| 
for a system of particles in quiescent fluid between the 
walls. The many-particle translational-rotational mobil- 
ity matrix fifj B is the inverse of the corresponding multi- 
particle resistance matrix Cil B ; i- e -> 



N r 

E 
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where I is the identity tensor. 

In our recent publications [2(| l27l |2^| we have intro- 
duced a formalism that allows us to efficiently evaluate 



A. Induced-force formulation 

In our approach, the effect of the suspended particles 
on the surrounding fluid is represented in terms of the 
induced-force distributions on the particle surfaces 



Fi(r) = a- 2 S{r 2 -a)t l {r) 



(8) 



where 



and r i = \r l 
field 



Ti=T- Hi (9) 

By definition of the induced force, the flow 



N 

v(r) = v oxt (r) + E / T ( r < r ') ' F *( r ') dr ' 

i—1 



(10) 



is identical to the velocity field in the presence of the 
particles [Mill HI- Here 

(11) 



T(r,r') = T (r-r') + T'(r,r') 



is the Green's function for the Stokes flow in the wall- 
bounded system, To(r) is the Oseen tensor (free-space 
Green's function), and T'(r,r') describes the flow re- 
flected from the walls. 

For a system of particles moving with the translational 
and angular velocities U, and !T; in the external flow 



the induced-force distribution 



can be obtained 



from the boundary-integral equation of the form 

AT 

[Zr 1 F,](r)+^ /[(l-^)To(r-r') 

j=i J 

+ T'(r,r')]-F J (r')dr'-vf(r)-v oxt (r), 



r G Si, (12) 
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where the rigid-body velocity field 

vf(r)=U 1 |ft,Xr 1 (13) 

and the external flow field v ext (r) are evaluated on the 
surface Si of particle i. In the boundary-integral equa- 
tion iJUJl, Zj denotes the one-particle scattering operator, 
which is defined by the relation 

F i = -Z i (v? l -vJ b ) ) (14) 

where vj n is the velocity incident to particle i. For specific 
particle models, explicit expressions for the operator Z^ 

are known [2i|HllH3- 

The force and torque acting on particle i can be eval- 
uated from the induced-force distribution using the inte- 
grals 

Ti = J F,(r)dr, T, = J r, X F,(r)dr. (15) 

The friction matrix © can be computed by solving the 
boundary equation (|12|l with the external flow v cxt in the 
form and no rigid-body motion, vj' b = 0. Similarly, 
the translational-rotational friction matrix is obtained 
by solving l|12|) with a nonzero rigid-body motion 113|) 
and no external flow, v cxt = 0. 

B. Multipolar expansion 

In our approach, the boundary-integral equation l|12|) 
is solved after transforming it into a linear matrix equa- 
tion. The transformation is achieved by projecting Ijl2|) 



The first term on the right side of the above expression 
corresponds to the single-particle scattering operator Z^ 
in equation |Q ; the second one to the integral operator 
with the kernel To, and the third one to the integral oper- 
ator with the kernel T'. Explicit expressions for the first 
two terms were derived by Felderhof and his collaborators 
[32l Isfj l3^ | some time ago. Quadrature relations j^E |2]J 
and asymptotic formulas |29| for the wall contribution 
Gjj were recently derived by our group (as discussed in 
Sec. HTTP] below). 



C. Friction and mobility of spheres in parabolic 
flow 

In order to evaluate the resistance tensors £* p and C 2 ip 
for immobile particles in Poiseuille flow, Eq. (|18fl is solved 



onto a spherical basis of Stokes flows. We use here the 
multipolar representation introduced by Cichocki et al. 
|32j | , but we apply a different normalization to emphasize 
full symmetry of the problem ptj [271 • 

Accordingly, the induced-force distributions at the sur- 
faces of particles i — 1 . . . N are expanded using the basis 
set of multipolar force distributions w^ n(T (ri). Similarly, 
the flows incoming to each particle are expanded into the 
nonsingular basis set of Stokes flows (rA Here I and 
m are the angular and azimuthal spherical-harmonics or- 
ders, and a = 0,1,2 characterizes the type of the flow. 
Explicit definitions of the basis sets wf ma and vf mcr (as 
well as their counterparts Wj~ i(T and v,~ _ that correspond 
to singular Stokes flows) are given in [2g, H2| • 

In order to obtain the multipolar representation of the 
boundary- integral equation (|12|l . we apply the multipolar 
expansion 

F((r) - ]T Mlm<7)a- 2 S( n - o)w+ a (r() (16) 

Ima 

to the induced- force density ©. The external flow rela- 
tive to the particle motion is similarly expanded, 

v rb (r) _ v e X t (r) = Ciilma^tmM). (17) 

Inserting these expansions into Eq. l|12f) yields a linear 
equation of the form 

N 

M i:j (lma | I'm'a^fjQ'm'a') = a(lma), (18) 

j — 1 I'm' a' 

where the matrix M can be decomposed as 



(19) 

I 

with the right-hand side corresponding to the velocity 
field Q). The resulting induced- force multipolar distri- 
butions l|16l) are projected onto the total force and torque 
using expressions 1)15(1 . The solution can be conveniently 
expressed in terms of the grand friction matrix 

F = M -1 , (20) 

which is inverse to the grand mobility matrix M with the 
elements given by Eq. (|19|) . 

As shown in |26j, the translational-rotational friction 
matrix (A, B = t, r) is given by the relation 

Cy B = X) X ( A I lm °) F v( lm ° I l'm'a')X(l'm'a' \ B). 

I'tna I'm' a' 

(21) 

Here F.ij(lma \ I'm' a') are the elements of the grand 
friction matrix (|2*U|) . and X(A | Ima) = X.*(lma \ A) are 



I 

M lj {lma | I'm! a') = SijSwS^Z^^a | a') + (1 - 8 l3 )G a lJ {lma \ I'm'a') + G'^lma \ I'm' a'). 
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the elements of projection matrices onto the subspace 
of translational (I = 1, a = 0) and rotational (/ = 1, 
a = 1) rigid-body motions. Explicit expressions for these 



matrices are listed in Appendix B of Ref . 129 . 

The resistance coefficients <^ Ap (A = t, r) are given by 
a relation analogous to (|2"T|l . 



J 



JY 



CfP = zZ X ( A I lm °) F ij(l™°- I I'm'a^Yjil'm'a' | p), 

7=1 Imcr I'm'cr' 



(22) 



r 



where Yj(l'm' a' \ p) are the elements of the matrix rep- 
resenting the orthogonal projection onto the subspace of 
pressure-driven flows J[J. Relation IL'l't and explicit ex- 
pressions for the matrix Yj(l'm'a' \ p) are derive in Ap- 
pendix E] We note that, unlike Eq. l(2"T]) . relation JUJ) 
involves summation over the particles. This summation 
is needed because the external parabolic flow is ap- 
plied to all particles in the system. 



D. Cartesian representation 

To determine the resistance coefficients l|21[) and (|22|l . 
the matrix l|19|l in the force-multipole equation l|18|l has 
to be first evaluated. Explicit expressions for the single- 
particle scattering matri x Z~ x and the free-space con- 



tribution G\j are known [32[ |3g • To evaluate the wall 
contribution GL to the matrix l|19(l we employ our re- 
cently developed Cartesian-representation method (2(|. 
For sufficiently large intcrparticle separations appropri- 
ate asymptotic expressions |2^ can also be used. 

As explained in pjj. [27J , the Cartesian-representation 
method relies on transformations between the spherical 
basis sets of Stokes flows v ; ^ ncr and the Cartesian basis 
sets v^ CT (where k is a lateral wave vector). According 
to the discussion in Sec. 1111 CI the multipolar spherical 
sets v^ ct correspond to an expansion of the velocity field 
into spherical harmonics. Due to symmetry, the matrix 
Zi, describing interaction of the flow field with a spherical 
particle, is thus diagonal in the spherical-harmonics or- 
ders I and to. The Cartesian basis sets correspond to an 
expansion of the velocity field into lateral Fourier modes. 
In the Cartesian representation the matrix Z w that de- 
scribes interaction of the flow with a wall is diagonal in 
the wave vector k. This diagonal structure of the scatter- 
ing matrices Zi and Z w yields a significant simplification 
of the problem. 

To express our results in a compact form, we intro- 
duce a matrix notation in the three-dimensional linear 
space with the components corresponding to the indices 
(7 = 0,1,2 that identify the tensorial character of the 
basis flow fields vrL . In this notation, a column vec- 
tor with components a(a) is denoted by a, and a ma- 
trix with elements A(a \ a') is denoted by A. Ac- 
cordingly, the column vectors associated with the coef- 



ficients fi(lma) and calmer) are represented by fj(Zm) 
and Ci(lm), and the two- wall Green's matrix with the ele- 
ments G'^ilma | I'm' a') is represented by G'y(Jm \ I'm'). 
We will also use 3x6,6x6 and 6x3 matrices composed 
of 3 x 3 blocks, as indicated below. 

Our result for the wall Green's matrix G', can be ex- 
pressed in terms of the following Fourier integral 



G'ijilm | I'm') = 



dk*(k;Z l ,Z„J?) 



(23) 



where g 



y — ^Hjcs, -r j,jcj is the projection of the vector 
R,j = Hi — Hj onto the x-y plane, and k = k x e x + k y e y 
is the corresponding two-dimensional wave vector. The 
matrix ^ in the above expression depends on the wall 
separation H and the vertical coordinates Zi and Zj of 
the points i and j (measured with respect to the position 
of the lower wall). This matrix is a product of several 
simple matrices, 

^(k-Z t ,Zj,H) = -7 7 - 1 T S c(;TO,k).S lW (k)-ZTw(k) 
• S Wj (k)-T cs (k,ZW), (24) 



where r\ is the fluid viscosity. 
The component matrices 



Tcs(Mm) = [Tsc(lm,k)]t = 



T+ s (k, Im) 
T cs"( k > lm ) 



(25) 



describe the transformations between the spherical (S) 
and Cartesian (C) basis fields. The superscripts ± refer 
to the singular and nonsingular basis fields for the spher- 
ical basis, and the fields that exponentially grow (+) or 
decay (— ) in the vertical direction z for the Cartesian 
basis. The matrices (|25|l consist of two 3x3 blocks cor- 
responding to the lower and the upper wall, respectively. 
Next, the matrices 



Sw s (k) = [S sW (k)] 



t _ 



'(kZ Ls ) 











~{kZ Vs ) 



(26) 

correspond to the propagation of the Cartesian flow-field 
components between the point s — i,j and the lower (L) 
and upper (U) walls. Here Z^ s = —Z s and Z\j s = H — Z s 
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E. Far-field asymptotic form 



Sc (R-iU ).••'' 



Sc(Ruj) 



Sc(R iL K 



f 

Ztw 



•••"Sc(R Lj ) 



FIG. 1: Schematic representation of Eq. 112 1 ti . The vectors 
Raj and RiA represent the relative position of the particle 
k = i,j and the lower (A = L) or upper (A = U) wall. 



are the relative vertical coordinates of the point s with 
respect to the walls. 
Finally, the matrix 



-TW 



(k) = 



Z- 1 S+ + (-kH) 



(kH) 



(27) 



describes scattering of the Cartesian flow components 
from the walls. The 3x3 matrices Z w represent scat- 
tering of the flow from a single wall, and the matrices 
S(t + (— kH) and S>Q~(kH) show the propagation of the 
flow field between the walls during the multiple-reflection 
process. The structure of the expressions (|23[) - (|27|1 is 
schematically represented in Fig. ^ The explicit expres- 
sions for the component matrices (|24l) - l|27() are listed in 
Appendix IT31 

We note that due to symmetries of the transformation 
and displacement matrices and the symmetry 



-TW 



(k) 



-TW 



(k)] t 



(28) 



of the two- wall scattering matrix, the Lorentz symmetry 
of the two- wall Green's matrix (|23|) is explicit. 



The exact Cartesian representation (I23[l -H27 [) of the 
wall contribution to the Green's matrix G', involves a 
two-dimensional Fourier integral, which has to be evalu- 
ated numerically. However, for sufficiently large interpar- 
ticle separations the calculation can be greatly simplified 
by using the far-field asymptotic expressions derived in 
|29j. Below we summarize this result. 

The derivation of the asymptotic expressions for the 
Green's matrix 

<^-=G°. + <4 (29) 
relies on the observation that for large lateral interpar- 
ticle distances, p\2 ^> H, the disturbance flow scattered 
from the particles assumes the Hele-Shaw (i.e., lubri- 
cation) form. Accordingly, the the far-field disturbance 
flow v as is driven by a two-dimensional harmonic pres- 
sure field p as , 



)V||j/ 



(30) 



The pressure p as is independent of the vertical variable z 
and satisfies the two-dimensional Laplace's equation 



Vjjp as (p) = 0, 



(31) 



where p — (x,y) is the lateral position with respect to 
the particle, and V|| is the two-dimensional gradient op- 
erator with respect to the lateral coordinates. The result 
(|30|l can be obtained using a lubrication expansion of 
the Stokes equations in the small parameter Hj p (where 
P=\p\)- 

To obtain the asymptotic expression for the Green's 
matrix G^ we use the results listed in Appendices lAl and 
IbI Accordingly, the asymptotic flow produced by a force 
multipole (|B3|> centered at the position of particle j is 
expressed in terms of the Hele-Shaw basis IjAljl using 
relation (|B2(1 . The resulting Hele-Shaw multipolar flow 
is translated to the position of particle i using the dis- 
placement formula (|A3|) . Finally the Hele-Shaw field is 
transformed back into the spherical basis using relation 
(|B1() . The above procedure [2!j yields a compact expres- 
sion of the form 



G a f (Zma | I'm' a') 



irr/H 3 



C(Zi;lma)S+Y{Qij;m \ m')C{Z 3 -l'm'a'), 
I 



(32) 



where the component matrices C and S^ yl are given by 
Eqs. jXPl and (|B4|l - l|B7|l . As explained in [^cj, the cor- 
rection 



5dj — dj — G as 



(33) 

to the asymptotic result l|32|) decays exponentially with 



the lateral interparticle distance g+j on the lengthscale 
H. Typically, the asymptotic approximation Gy ~ G as 
yields accurate results for Qij/H > 2. 
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d/H_ 

Z/a 0.999 0.995 0.990 0.950 0.900 0.500 0.200 

1.1 0.641 0.583 0.286 

1.01 0.418 0.498 0.520 0.401 0.188 

1.007 0.415 0.483 0.502 0.382 0.179 

1.005 0.376 0.409 0.469 0.486 0.366 0.171 

1.001 0.304 0.350 0.368 0.409 0.419 0.306 0.141 

TABLE I: Normalized translational velocity U x /U p of a single sphere of diameter d = 2a in imposed parabolic flow for 
different wall separations H and particle positions Z with respect to the lower wall. 







d/H 










Z/a 


0.999 0.995 


0.990 


0.950 


0.900 


0.500 


0.200 


1.1 








0.0197 


0.723 


1.189 


1.01 




5.14E-4 


0.101 


0.177 


0.620 


0.903 


1.007 




0.0159 


0.109 


0.181 


0.600 


0.866 


1.005 


1.95E-4 


0.0269 


0.115 


0.184 


0.582 


0.834 


1.001 


2.34E-5 0.0362 


0.0556 


0.127 


0.183 


0.504 


0.705 



TABLE II: Normalized angular velocity HQ y /U p of a single sphere of diameter d — 2a in imposed parabolic flow Q, for 
different wall separations H and particle positions Z with respect to the lower wall. 



F. Numerical implementation 

In order to determine the resistance matrices <|21[) and 
1)220. the induced- force- multipole equation l|18l) is solved 
with the matrix (|19f) evaluated using known results |32l 
l38l | for Z~ x and G?-, and our Cartesian representation 
(|23|l for G\i . For sufficiently large interparticle distances 
a simpler relation i|32|) may be be used instead. After the 
friction matrices have been obtained, the mobility matrix 
(0 can be calculated from expressions © and J7J. 

To accelerate numerical convergence of the Fourier in- 
tegral in (|23|) (especially, when both particles i and j are 
close to a single wall), the integrand (|2~4"|) is decomposed 
into two single-wall contributions and $u and the 
correction term 

*(ft)=¥ L (A) + *u(fc) +<**(*)• (34) 

The single wall contributions can be integrated analyti- 
cally [39L |4C| . Moreover, as shown in [2|j , the correction 
term 5^(k) is easier to integrate numerically than the 
original highly oscillatory integrand ^(k). 

As in other numerical algorithms based on a multipo- 
lar expansion of Stokes flow [3^, the force-multipole 
equation QlSj l is truncated at a given multipolar order 
I = 'inns before it is solved numerically. To accelerate the 
convergence of the approximation with i max we employ 
the standard lubrication correction |42( on the friction- 
matrix level. We closely follow the implementation of 
the method described in [3^ (for a single wall problem). 
Accordingly, the translational-rotational friction matrix 
Cij = Cij B {A, B = t, r) is represented as a combination 

C, C/- 2 • • AC, (35) 



of the two-particle superposition contribution in free 
space d S j UP ' 2 j the single-particle/single- wall superposition 
contribution £ f . w , and the remainder A^. The super- 
position contributions Ci/ P ' 2 an d CiJ P,W are determined 
very accurately using the power-series expansions of the 
friction matrix in the inverse interparticle separation and 
the inverse distance between the particle and wall, re- 
spectively. The remainder A£jj , evaluated as a difference 
between the multipolar expansion of the full friction ma- 
trix and the superposition contributions, converges with 
^max much faster than the full friction matrix Cj itself. 

In the present implementation of our method, the lin- 
ear equation (|19|) is solved by matrix inversion. Thus, 
the numerical cost of the calculation scales as 0(N 3 ) 
with the number of particles N. (Numerical cost of this 
order is typical of unaccelerated Stokesian-dynamics al- 
gorithms.) We note, however, that the PPPM or fast- 
multipole acceleration techniques |43l ] can naturally be 
used in our Cartesian-representation algorithm — we will 
return to this problem in our future publications. 

IV. RESULTS AND DISCUSSIONS 

We now present some characteristic examples of single- 
and many-particle results. We consider both the mo- 
tion of freely suspended particles in the external flow 
and forces and torques on fixed particles subjected 
to this flow. The results for an isolated particle are 
obtained with the truncation at the multipolar order 
Imax = 32, which yields accuracy better than 0.1 %. For 
two-particle and multi-particle systems we use l ma x = 12 
and Z max = 8, respectively. The corresponding accuracy 
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is of the order of 1 %. 



A. Single particle system 

Motion of a single particle in a parabolic flow between 
two planar walls was recently considered by Staben et 
al. |13| and by Jones (2^| (see also much earlier results 
by Ganatos et al. 01)- We thus give here only limited 
results for this system. In Tables [I] and [H] we list a set 
of our highly accurate results for the linear and angular 
velocities (JSJ of a force- and torque-free particle in the 
parabolic flow The linear velocity U is normalized 
by the magnitude of the parabolic flow U p , and the an- 
gular velocity l~i by U p /H. Only the x component of the 
linear velocity and the y component of the angular ve- 
locity are given because all the other components vanish 
by symmetry. 

In order to verify our results and test the accuracy of 
the calculations reported in [^3], the velocities U x and 
VL V are given for a subset of configurations represented in 
Tables I and II of [30J. We also present some additional 
results for tight configurations with H w d. 

We find that our results are in good agreement (up to 
three digits) with those reported in [3J| for Z/a> 1.01, 
where Z is the position of the particle center measured 
from the lower wall. For smaller gaps between the wall 
and the particle the discrepancies are about 1.5%. An ex- 
ception is the rotational velocity in the tightest configura- 
tion reported in H3 (i.e., d/H = 0.95 and Z/a = 1.007), 
where the error is 11%. We expect that these discrep- 
ancies stem from inaccuracies of the boundary-integral 
calculations in [30| — the convergence tests we have per- 
formed indicate that the accuracy of our results is bet- 
ter than 0.05%. We also note that our results agree 
with those of Jones [2^| and with our earlier results of 
a multiple-reflection method |45| . 



B. Two- particle system 

1. Particle velocities 

A sample of characteristic results for the translational 
and rotational velocities Uj and fi, (i=l,2) of two force- 
and torque-free particles in the parabolic flow Q are pre- 
sented in figures [21[SJ The linear and angular velocities 
are normalized in these plots by U p and U p /d, respec- 
tively. The results are plotted versus the lateral particle 
distance pyijH for a moderate channel width 

H = 2d. (36) 

In Figs. |3 and 0] the particle 1 is in the center position 
Z\ = H/2, and in Fig. [3] and [S] it is in the off-center posi- 
tion Z\ = H/3. The results are given for several vertical 
positions Zi of particle 2. (We recall that Zi denotes the 
distance of particle (i) from the lower wall.) In Figs. [21 



and |3 the particle pair is oriented in the longitudinal di- 
rection x and in Figs.0]and[S]in the transverse direction y 
with respect to the flow. We note that U y = fl x = fi* = 
for the longitudinal configuration and U y = U z = VL X = 
for the transverse configuration, by symmetry. 

The results in Figs. indicate that the effect of mu- 
tual particle interactions is small if both particles are at 
the same vertical position in the channel. The effect is 
the largest if one of the particles is near the channel cen- 
ter and the other close to a wall. The results also reveal 
a different behavior in the near-field and far field regions, 
as discussed below. 

a. Near- contact and intermediate region The results 
in Figs. indicate that the dependence of the linear 
and rotational particle velocities on the interparticle dis- 
tance is much more complicated in the wall-bounded sys- 
tem than in free space. This complex behavior stems 
from the competition between the tangential and normal 
lubrication forces and backflow effects associated with 
the velocity field scattered from the walls. 

For near-contact particle configurations 

£i2 « 1 (37) 

(where ei2 = Ri^/d— 1 is the dimensionless gap between 
the particle surfaces, and R\2 = |Ri — R2I) the particle 
dynamics is strongly influenced by the lubrication forces. 
The normal relative particle motion is arrested by the 
0(612) normal lubrication force at the dimensionless gap 
612 of several percent. The relative tangential and rolling 
motions are opposed by much weaker O (log £12) lubrica- 
tion forces. These motions are thus still quite substantial 
for ei2 ~ 10 -3 and vanish only for nonphysically small 
gaps. 

A decrease in the relative tangential and rotational 
particle motion at small interparticle distances results 
in an increased overall dissipation, which may cause a 
decrease of the horizontal particle velocities even in sym- 
metric particle configurations with Z2 = Z± or Z2 = 
H-Zi (cf., the results for Z 1 /H = § and Z 2 /H = |, § 
111 Fig. EJ). We note that a pair of touching particles in a 
transverse configuration (Figs. 01 and 0) does not move, 
in general, as a rigid body, because there is no lubrication 
resistance to the relative particle rotation around the axis 
connecting their centers. 

In some cases the normal and tangential lubrication 
forces have an opposite effect on a given velocity com- 
ponent. This produces sharp kinks in some curves at 
near-contact positions (e.g., U\ z and U% x for Z 2 = 2.67a 
in Fig. EJ). An additional change of sign of particle ve- 
locities relative to the velocities at infinite interparticle 
separations p\2 — > 00 may occur due to a backflow as- 
sociated with scattering of the flow from the walls. Due 
to a combination of the lubrication and back-flow effects, 
the z component of the particle velocities changes sign 
twice for some longitudinal configurations. 

b. Far-field region As discussed in Sec. IIII El for 
large lateral interparticle distances, the hydrodynamic 
interactions in a wall-bounded system are determined 
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FIG. 2: Normalized linear and angular velocities of two spheres of diameter d = 2a in imposed parabolic flow versus 
particle lateral distance pi2 normalized by wall separation H . The particle pair has the longitudinal orientation with respect 
to the flow direction, the wall separation is H — 2d, and sphere 1 is in the center position Zi = 2a. The positions of sphere 2 
are Z2 = 1.02a (short-dashed line), Z2 = 1.33a (long-dashed line), Z2 = 2.0a (solid line). 



by the far-field form of the disturbance flow scattered 
from the particles. The scattered flow has the Hele-Shaw 
form described by Eqs. (|3(J|I and (|31f) . We recall that the 
asymptotic form of the flow field is approached exponen- 
tially on the lengthscale H . 



Taking into account the symmetry of the problem we 
find that the far-field disturbance velocity produced by 
a particle in external flow Q is given by equation i|3U|) 



with the pressure of the form 



cos< 



P 



(38) 



where <f> is the polar angle between the lateral position 
vector p and the flow direction e x . To the leading order in 
the multiple scattering expansion, relations l|30|l and 138|) 
determine thus the far-field form of the hydrodynamic 
resistance and mobility functions for a pair of particles 
in the external parabolic flow. 
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FIG. 3: Same as figure [2] except that particle 1 is in an off-center position Z\ = 1.33a, and the positions of particle 2 are 
Z 2 = 1.02a (short-dashed line), Z 2 = 1.33a (long-dashed line), Z 2 = 2.0a (solid line), Z 2 = 2.67a (dotted line), Z 2 = 2.98a 
(dash-dot line). 



In particular, Eqs. (|3U[) and (|38[) indicate that the flow 
v as has only lateral components. It follows that the z 
components of the translational and rotational particle 
velocities (@J vanish in the far-field regime. This behav- 
ior is clearly seen in Figs.|2HSl where these velocity com- 
ponents approach zero exponentially. 

Next, the disturbance field (|3U[1 with the pressure given 
by Eq. (|3"5|) behaves as 

v as ~ \. (39) 
P 

Thus the linear and angular lateral velocities shown in 



Figs. approach the one-particle asymptotic values as 
0(p~ 2 ). The result Ip^ should be contrasted with the 
behavior 

v as ~ - (40) 
P 

in free space and 

v as ~ \ (41) 
P A 

in the presence of a single wall (where we assume that 
Z\,Zi < p). According to Eqs. P? )h |2T )1 . the decay 
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FIG. 4: Same as figure [2] except that for the transverse orientation of the particle pair. 



of the far-field flow in the presence of one or two walls is 
faster than the corresponding decay in free space, because 
the walls absorb momentum, and thus they slow the fluid 
down. 

On the other hand, the decay of the flow field (|39|l in 
the two-wall system is slower than the decay (|41fl in the 
presence of a single wall. This behavior stems from fluid- 
volume conservation constraint. In the system confined 
by a single wall the fluid displaced by the particle pri- 
marily flows above the particle, far from the wall, where 
it encounters small resistance. In contrast, in the pres- 
ence of two walls, the flow is limited to the quasi-two- 
dimensional region; hence, it has a longer range. 



Since the total flux associated with the quasi-two- 
dimensional flow l|39fl vanishes for p — > oo, the fluid veloc- 
ity must form a backflow pattern, unlike the behavior in 
the unbounded three-dimensional space. The backflow, 
described by the dipolar velocity field lj3UI) and re- 
sults in an enhancement of relative particle motion for 
the transverse orientation of the particle pair (as seen for 
some particle configurations in the top panels of Figs. 01 
and [SJ . We note that an analogous behavior was dis- 
cussed by Cui et al. [19| in their study of pair diffusion 
in a confined, quasi-two-dimensional colloidal suspension. 
Similar effect was also independently described in our re- 
cent papers j2(| |2tJ . 
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FIG. 5: Same as figure |^1 except that for the transverse orientation of the particle pair. 
2. Crossover behavior force perturbation 



The far-field disturbance flow, discussed above, affects 
not only the velocities of freely suspended particles, but 
also the hydrodynamic resistance force © acting on im- 
mobile particles in the external flow (JIJ. Figures [|J] and 
illustrate the crossover of the resistance force between 
the three regimes corresponding to the disturbance flows 
of the form given by equations i|39|l - H41[) . To emphasize 
the behavior of the force in the far-field regime, the re- 
sults are shown for the x component, SF^, of the rescaled 



where 



SFi 



{Ti T?)/F st , 



(42) 



(43) 



with F st — 37T7/G? denoting the Stokes resistance force, and 
ZF?° representing the value of the force J~i for pi2 — > 00. 

In Figs. El and force perturbation SF ix is plotted ver- 
sus the lateral particle separation P12 for two particles at 
the same vertical position Z\ — . In one configuration, 
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FIG. 6: Rescaled force perturbation (1421 versus interpar- 
ticle distance pi2 normalized by the particle diameter d for 
a pair of particles in longitudinal orientation. Wall separa- 
tion H = 1.02d (solid line), H = 5d (long-dashed lines), 
H — 16d (short-dashed lines), H — oo (dotted lines). The 
top three lines correspond to particles in the center plane 
Z\ — Z2 = \H, and the bottom three to particles in the near- 
wall configuration Z\ — Z2 = 1.02a. For H — 1.02d (middle 
line) the center and near-wall configurations coincide. 

the particles are at the center plane 

Z 1 = Z 2 = \H, (44) 

and in the other one they are close to the lower wall, 

Z X = Z 2 = 1.02a. (45) 

The results are shown for several different wall separa- 
tions. Since the particles are at the same vertical posi- 
tion, the force 5F X = SFi X is independent of the particle 
index i. Figure |H| represents the results for the longitudi- 
nal orientation of the particle pair, pi 2 — pi 2 e x , and Fig. 
0the results for the transverse orientation p\ 2 = pi 2 e y . 
The force perturbation (|42|l for the longitudinal orien- 
tation is positive. It is shown on the logarithmic scale 
to emphasize the algebraic asymptotic behavior. For the 
transverse orientation the perturbation force in the wall- 
bounded systems changes sign due to the backflow effects 
discussed in Sec. IIVB II The results are thus plotted on 
a linear scale in two separate panels for the center (top 
panel) and the near-wall (bottom panel) configurations. 

The results shown in Figs. El and {7\ clearly demonstrate 
the crossover between different regimes corresponding to 
the far-field disturbance velocity fields of the form 
l)4ip. For very large wall separations H — ► 00 and the 
center particle position (|!*4f> the rescaled force pertur- 
bation H42|) behaves as SF X ~ pi%, which indicates that 
5 F x = 0(pi~2 ) i m agreement with the estimate <|4(J|) of 
the disturbance-flow magnitude in free space. For the 
near-wall position l|45|l and the longitudinal orientation 
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FIG. 7: Same as Fig. |H1 except that for the transverse orienta- 
tion of the particle pair. The top panel corresponds to parti- 
cles in the center plane Z\ — Z2 = \H, and the bottom panel 
to particles in the near-wall configuration Z\ = Z2 — 1.02a. 

of the particle pair we find SF X = 0(pi 2 3 ), consistently 
with the estimate Q4ip. In contrast, 8F X = 0(p^ 2 5 ) m the 
transverse case, due to an additional cancellation of the 
far- field contributions. 

For finite wall separations the force perturbation 
crosses over from the above-described behavior in the 
regime a <C P12 <C -ff to the behavior 5F X = 0(pi 2 2 ) 
(i.e., SF X ~ const) for pi 2 3> H, in agreement with Eq. 
(|39|l . Typically, the far-field behavior 0(pi 2 ) is observed 
already for p\ 2 > 2H. 



C. Multiparticle systems 

In this section we examine the influence of the walls on 
the hydrodynamic interactions in confined multi-particle 
systems. We focus on collective phenomena that involve 
cumulative effects of the far-field flow l|39|l . As shown in 
our recent studies of particle motion in quiescent fluid 
[26l I27I 129) . the backflow associated with the dipolar 
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form l|30(l and l|38(l of the far-field velocity may produce a 
strong, positive feedback resulting in large magnitudes of 
induced forces. In such cases the far-field flow dominates 
the behavior of the system. Below we examine similar 
phenomena for particles in the imposed parabolic flow. 



1. Motion of linear arrays of spheres 

First we analyze the effect of confinement on the mo- 
tion of rigid linear arrays of touching spheres. In earlier 
papers j2y, [2?], H^J , we have shown that the behavior of 
such arrays in quiescent fluid is strongly affected by the 
walls. In particular we have demonstrated that, unlike 
in free space, the hydrodynamic resistance force in chan- 
nels with H sa d depends significantly on the orientation 
of the array with respect to its velocity. If the orienta- 
tion of the array, moving along the channel, is parallel to 
the velocity, the resistance force evaluated per one sphere 
decreases with the length of the array. In contrast, for 
the transverse orientation the resistance force per particle 
increases nearly linearly with the array length. This in- 
crease results from the pressure buildup associated with 
the positive-feedback backflow effects. 

The motion of linear arrays of spheres in the imposed 
parabolic flow (JTJ is illustrated in Fig. [S] and [5J The 
arrays are parallel to the walls and are oriented either in 
the longitudinal direction x or the transverse direction y. 
Figure|Hlpresents the translational velocity of arrays with 
different length, placed in the mid-plane z = \H. The 
results are given for several channel widths H. Figure 
shows linear and angular velocities of arrays at different 
vertical positions in the channel. The linear velocities are 
non-dimensionalized by the local velocity of the imposed 
flow 



v 



w cxt (Z) 



and the angular velocities by the local share rate 



7o = 



dv cxt (Z) 
dZ 



(46) 



(47) 



evaluated at the position Z of the array center. For the 
mid-plane position z — ^H, vq is identical to the ampli- 
tude Up of the imposed flow QJ. 

The results in Fig.[S]indicate that the normalized veloc- 
ity of an array U x /U p is smaller in channels with smaller 
width. This behavior stems primarily from the increased 
dissipation in the gaps between the particles and the 
channel walls. The decrease of the mobility is strongest 
for long arrays in longitudinal orientation — the far-field 
disturbance flow produced by each of the particles op- 
poses the motion of the array in this case. For the trans- 
verse orientation, the scattered flow acts in the direction 
of the external flow; due to the cooperative feedback ef- 
fects longer arrays move faster than the shorter ones. In 
narrow channels with H sa d, very long chains in trans- 
verse orientation translate with the velocity that is close 
to the average velocity of the unperturbed fluid. 



A set of results for short (N = 3) and a long (N = 20) 
linear arrays at off-center positions in channels with dif- 
ferent width is presented in Fig.(SJ The configurations are 
parametrized by the normalized distances of the particle 
surfaces from the lower and upper walls, 



= \{Z-a)/a, e v = \(H-Z-a)/a. (48) 



The translational and rotational velocities are shown for 
arrays at two vertical positions £l = 0.001 and 6l = 1.1, 
and they are plotted versus the distance eu- 

The results in the upper panels of Fig. indicate that 
the translational velocity of an array at a fixed distance 
from the lower wall diminishes rapidly with the decreas- 
ing eu <C £l due to the O(logeu) lubrication resistance 
associated with the interaction with the upper wall. In 
the case of the longitudinal orientation of the chain, the 
translational and rotational velocities saturate at eu ~ 1 • 
In contrast, for the transverse orientation, the effect of 
the upper wall on the translational velocity of the ar- 
ray has a much longer range, especially for the larger 
value of the chain length N. Moreover, the effect of the 
upper wall is more pronounced for el = 0.001 than for 
e = 0.1. These observations are consistent with the back- 
flow mechanisms discussed above. 

Lower panels of Fig.|5|represent the normalized angular 
velocity Oy/70 of the arrays. We note that the angular 
velocity itself changes sign for = e\j; however, the nor- 
malized quantities shown in Fig. [5] are positive, due to 
the corresponding change of sign of the local share rate 
(|47|l . For the longitudinal orientation of the chain the 
angular velocity is several orders of magnitude smaller 
than the angular velocity in the transverse case. This 
strong effect can easily be explained in terms of particle- 
wall lubrication forces. The rotation of the chain oriented 
perpendicularly to the flow is governed by the O(logeL) 
lubrication resistance. The rotation of the chain oriented 
parallel to the flow involves motion of individual particles 
towards the wall and away from it, and thus the lubrica- 
tion forces are much stronger 0{e : L 1 ). 



2. Hydrodynamic drag on adsorbed particles 

In Figs. I1UH12I the results are presented for the hydro- 
dynamic friction forces and torques on individual parti- 
cles in arrays of spheres adsorbed on one of the walls in a 
parallel- wall channel. Understanding of such forces is im- 
portant in an analysis of the removal of colloidal particles 
from a wall by an applied flow 0]. Moreover, our results 
provide a further illustration of hydrodynamic phenom- 
ena associated with the far-field form of the disturbance 
velocity field produced by the particles. 

Figures EH an d IHI show forces and torques acting on 
individual particles in closely packed arrays of touching 
spheres. Arrays that are loosely packed are considered 
in Fig. 1121 The results are given for a single sphere, lin- 
ear arrays of spheres, and hexagonal arrays of spheres. 
The horizontal components of the forces and torques are 
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FIG. 8: Normalized translational velocity U x /U p of rigid, linear arras of touching spheres, for the longitudinal (left panel) and 
transverse (right panel) orientation with respect to the imposed flow Q. The results are shown versus the number of particles 
in the array N. The arrays are in the center plane z = between walls separated by H = 1.004d (solid circles), H = 1.02d 
(open circles), H = l.ld (solid squares), H = 2.0d (open squares), H — 4.5d (crosses). 



represented by line segments and the normal components 
by the dashed (orientation away from the wall) and dot- 
ted (towards the wall) circles. The forces are scaled by 
arjvo and the torques by a 2 r]Vo, where vq is the local fluid 
velocity l|46|) . The lengths of the line segments and the 
radii of the circles are proportional to the magnitude of 
the represented quantity. 

The results presented in Fig. ^] indicate that the lat- 
eral drag force on a single particle only weakly depends 
on the channel width, while the forces on particles in lin- 
ear arrays vary almost by a factor of five when the wall 
separation H changes from 2.02a to infinity. Similarly 
strong dependence of the forces on the channel width is 
observed for the two-dimensional arrays of spheres. 

The hydrodynamic drag forces are the largest for lin- 
ear arrays of spheres in narrow channels of the width 
only slightly bigger than the particle diameter (cf., the 
top panel of Fig. fTOf) . The large forces are associated 
with the pressure buildup in front of the array (29j | . As 
explained in Sec. IIV C II the pressure buildup involves 
interaction of the flow with both walls in an essentially 
non-additive manner. Thus, as shown in [2fij . this effect 
is not captured by the usual approximation based on a 
superposition of two single- wall contributions. 

The results for two-dimensional arrays of spheres in- 
volve significant screening effects resulting from mutual 
particle interactions. Accordingly, the lateral forces act- 
ing on individual spheres in two-dimensional arrays are 
smaller than in the corresponding linear- array systems. 
The forces on the first and last row of particles are larger 
than the forces on particles near the center of the array. 
This effect is most pronounced for large wall separations 
H — in narrow channels the relative force differences are 
smaller due to the quasi-two-dimensional character of the 
flow. 

Our results for the vertical forces indicate that their 
magnitude is much smaller than the magnitude of the 



lateral forces. Indeed, the only significant vertical forces 
exist on the first and last rows of particles in two- 
dimensional arrays. The maximal value of these forces 
occurs for H as 4a. When the channel width is smaller, 
the upper wall suppresses the vertical flow. On the other 
hand, when the gap between the top wall and the spheres 
is too large, the volume of fluid deflected by the array is 
distributed over the larger space, and the vertical flow be- 
comes weaker. We note that there are no vertical forces 
on particles in linear arrays because of the flow-reversal 
symmetry. 

The behavior of the torque exerted on the particles by 
the fluid is illustrated in Fig. ^] Characteristic features 
of the torque distribution can be explained using argu- 
ments similar to the ones given above. For example, the 
lateral torques are the largest for H as 4a, for the same 
reason as the corresponding behavior of the normal force. 
Our results indicate that the vertical component of the 
torque is significant only for particles at the edges of the 
arrays, especially for small values of H . 

Figure^Jshows plots of forces on the spheres in loosely 
packed arrays of spheres. For linear arrays, the lat- 
eral forces are relatively small even for H as 2a because 
the flow can pass through the inter-particle gaps with- 
out building up a substantial pressure drop. Moreover, 
the forces on particles in different positions in such ar- 
rays are of approximately equal magnitude. For the two- 
dimensional loosely packed arrays the lateral forces are 
larger than for the closely-packed case, which indicates 
that the screening effects are smaller. 



V. CONCLUSIONS 

We have used our recent Cartesian-representation al- 
gorithm to study hydrodynamic interactions of spherical 
particles in a parabolic flow between two parallel planar 
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FIG. 9: Normalized translational and angular velocities of rigid, linear arras of touching spheres with length TV = 3 (dashed 
lines) and N = 20 (solid lines), for the longitudinal (left panel) and transverse (right panel) orientation with respect to the 
imposed flow Q. The distance 1481 of the particle surfaces from the lower wall eL, as labeled; the results are shown versus 
the distance from the upper wall etr. The translational velocity U x is normalized by the local fluid velocity 14611 : the angular 
velocity Q y is normalized by the local shear rate 14711 . For the longitudinal configuration the angular velocity is in additional 
rescaled by TV 2 . The inset shows the region with small e\j. 



walls. An important feature of our method is that at 
each multipolar-approximation level the boundary con- 
ditions at the walls are exactly satisfied. This ensures 
that the far-field flow produced by the particles has a 
correct form of a two-dimensional Hele-Shaw lubrication 
velocity field. Our analysis indicates that the far-filed 
flow and the associated backflow effects strongly affect 
hydrodynamic interactions in confined multiparticle sys- 
tems. 

We have presented a set of numerical results for a sin- 
gle particle, a pair of particles, and arrays of many parti- 
cles. Our o ne-p article calculations agree well with earlier 
results [2^, |3(J, . For very tight configurations with 
small gaps between the particle surface and the walls we 
provide more accurate data than those reported previ- 
ously poj. For two-particle and multi-particle systems 
no accurate results have been available so far. 

Our numerical calculations reveal that the pair and 
multiparticle hydrodynamic interactions in the wall 
bounded system are much more complex than the in- 
teractions in free space. In particular, unlike in free 
space, the sign of mutual friction and mobility functions 



depends on the relative particle position in the flow- 
vorticity plane. The changes of sign result form com- 
bined effects of the short-range dissipation in the near- 
contact regions and backflow due to the confinement. 
Related backflow phenomena were recently observed in 
quasi-two-dimensional suspensions of Brownian particles 
|l9j . For elongated particles in narrow slit pores, such 
a backflow results in a strongly non-isotropic diffusion 
constant pi WA. EME^ . 



The far-field flow also determines the fluid velocity 
distribution and the hydrodynamic drag forces in two- 
dimensional arrays of particles adsorbed on a wall. Our 
results indicate that in narrow channels with the width 
H similar to the particle diameter d the hydrodynamic 
forces act mostly in the horizontal direction. Normal 
forces, which may lead to particle resuspension, are max- 
imal in channels with H sa 2d. Our results on hydrody- 
namic drag on immobile absorbed particles can be used 
in an analysis of particle removal from a channel by a 
flow. 
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FIG. 10: Hydrodynamic drag forces on a single sphere and on 
individual particles in linear and hexagonal arrays of spheres 
adsorbed on a wall in a channel with the width H (as labeled). 
The spheres are depicted by solid circles. The lateral forces 
are represented by the line segments, and the normal forces 
by dashed (force away from the wall) or dotted (towards the 
wall) circles. A line segment (circle) of the length (radius) 
equal to the particle radius a represents a force of magnitude 
207ra?7?;o, where vq is the local fluid velocity 14611 . 
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FIG. 11: Same as Fig. 111)1 except that the results are for the 
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APPENDIX A: HELE-SHAW BASIS 

As shown in [23 , the far filed-flow in the two- wall ge- 
ometry has the Hele-Shaw, lubrication form. Such a flow 
can be represented in terms of singular (— ) and non- 
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H = l.Old 



o- or 



(|A2b[) and (|A1|) . The prime at the summation sign is 
introduced to emphasize that this term is omitted. 
The matrix 



S*, (e;m | ro') = 6(-mm>)(-ir i^,^,,, ' 



(A4) 

in Eq. I|A3(I is identical to the displacement matrix for 
the two-dimensional harmonic potentials JX2j. We note 
that due to the presence of the Heaviside step function 



6{x) = 



0, 
1. 



x < 0, 
x > 0, 



(A5) 



H = oo 
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FIG. 12: Same as Fig. 1101 except that the results are for 
loosely packed arrays of spheres. 



singular (+) Hele-Shaw basis fields of the form 



where 

®o(p) = -Inp, 
are singular and 



m[P> 2\m\ P 



^(p)=P |m| e 



(Al) 

(A2a) 
(A2b) 



non-singular two-dimension harmonic functions. Here p 
is the lateral position vector with the polar coordinates 
(p,4>), and Vii is the two-dimensional gradient operator 
with respect to the lateral coordinates. 

The Hele-Shaw flow fields IjAlfl centered at lateral po- 
sitions Qi and Qj are linked by the displacement formula 



<T(r - Qi) = J2' v « + ( r - Qi)S? y T(eif,™ I "0. 

m— — oo 

(A3) 

where Qij = Qi — Qj. The term with m = in the above 
relation vanishes because Vq S+ = according to Eqs. 



the Hele-Shaw basis fields with the same sigh of indices 
m,m' 7^ do not couple in the displacement relation 
31- 



APPENDIX B: TRANSFORMATION BETWEEN 
THE HELE-SHAW AND SPHERICAL BASIS 
SETS 

In this Appendix we list some relations between the 
Hele-Shaw basis of asymptotic far field flows (IA1|) and 
the multipolar spherical basis fields defined in 32] (in the 
normalization introduced in (2(|). 

As shown in j2(J , the nonsingular Hele-Shaw field + 
centered at the lateral position Qi has the following ex- 
pansion in terms of non-singular spherical basis fields 
v zm<r centered at R l = Qi + ZS Z , 

v- + (r - Qi) = £v+ m(T (r - R,)C(Z i; Una). (Bl) 



There is also a reciprocal expression 

QjiZj) = ~ Qj)C(Zj;lma) 

(B2) 

for the far-field flow ufL a produced, between the walls, 
by a force multipole 



nr) = a- 2 S(r j -a)w+ ma (r j ) 



(B3) 



(cf., the multipolar expansion (|16|> ^ . 

Explicit expressions for the transformation matrix 
C{Zi]lma) have been derived in [29j. Accordingly, the 
nonzero elements of C(Zf, Ima) satisfy the condition 



I + a - \m\ < 2, 
and they can be written in the form 



(B4) 



C(Z;l±p*) = B± lJia (Z;p J ), p=\m\> 1. (B5) 
Here, the B^ a (Z; p) are the elements of the 3x3 matrix 
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-Z{H - Z) 

-H(H-2Z) 
(^ + l)(2/i + 3)V2 

2/i(M + l) 1/2 
. (m + 2)(2^ + 3)(^ + 5) 1 / 2 

I 



T{H-2Z) 

±2/X 

(M + 1)(2 M + 3)V2 



(B6) 



with 



A^) = ( T 2)V 



4tt 



(2 M + 1)(2 M )! 



1/2 



(B7) 



The range A = 0, 1, 2 of the index X — I ~ \m\ in equation 
(IB6II result from the conditions Iml < I and 1B4I). 



APPENDIX C: PROJECTION MATRIX Y 

Relation l|22[l is derived by inserting expansion Q17I I of 
the external flow (HJ) into Eqs. (141) and (145) in Ref. |26|. 
Comparing Eqs. and (|22() to the resulting formula we 
conclude that the the matrix Yj(lma \ p) is determined 
by the expansion 

v cxt (r) = £ v+ mff (r - R,)r,(Wb) (CI) 



of the parabolic flow into the spherical basis centered 
at Rj. 

To obtain the explicit expression for the matrix 
Yj{lma | p), we represent the external parabolic flow 
in terms of the Hele-Shaw asymptotic basis (IA1I) . 

v cxt = _4 C / p JJ-2( v as+ + v aa+)_ ( C2 ) 

Inserting expansion (|Blfl into the above relation and com- 
paring the result to (|C1(I yields 

Yj(lma\p) = ~4U p H- 2 (S ml + 5 m - 1 )C(Z J ;lma), (C3) 



where <5 m fc denotes the Kronecker delta. Using relations 
(IB4ll-(IB7ll for the matrix C we thus find 



{Y 3 (\ + 1 ± 1 <7|p)} Aiff=0il , = -4H- 2 U f 



-ZjiH-Zj) T{H-2Z j ) 2 



-{H-2Z, 
2^5 

2 

15%/3 



±1 

7E 





(C4) 



All other elements of Yj vanish, by Eq. (|B4|I . [Sukalyan : 
verify !?.] 



APPENDIX D: COMPONENT MATRICES IN 
CARTESIAN REPRESENTATION 



In this Appendix we list explicit expressions for the 
component matrices in the Cartesian representation 
(I23|) 127|) of the Green's matrix G'ij. These expressions 
are derived in Ref. 26]. 

The component transformation matrices in Eq. I|25|) 



can be represented in the factorized form 

T± s -(k,H = [T+±(Zm,k)]t = r(27rk)-y 2 e im ^ 

xf± s -(Zm)-K(fc,0, (Dl) 

where (k, ip) are the polar coordinates of the vector k. In 
the above expression 



K{k, l; a | a ) = 5 aa 'k 



l+a-l 



and 



T 



cs 



(-1) 



c -2b 4a 
6 -2a 
a 



(D2) 



(D3a) 
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a 
b 2a 
c 2b 4a 



(D3b) 



The three independent scalar coefficients in equations 
(TED are 



[4(/-ra)!(Z + m)!(2Z + l)] 



-1/2 



(D4a) 



6 = 2am/l, (D4b) 
Z(2Z 2 — 2Z — 1) — 2m 2 (Z - 2) 



c = a- 



1(21 - 1) 



(D4c) 



The component displacement matrices in Eqs. <|26|) and 
(|2"7J) are given by the relation 



S++(kZ) = [S c -(-fcZ)]t = 



1 2kZ 
10 
1 



JtZ 



(D5) 



For rigid walls with no-slip boundary conditions the 
single-wall reflection matrix in Eq. (|27Jl has the form 



Z w = 



1 
1 
1 



(D6) 



For planar interfaces with other boundary conditions 
(e. g., a surfactant-covered fluid-fluid interface discussed 
in |48|) the scattering matrix is different from identity, 
and it may depend on the magnitude of the wave vector 
k. Explicit expressions for scattering matrices for such 
systems will be presented elsewhere. 
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